Fluctuations and oscillations in a simple epidemic model 
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We show that the simplest stochastic epidemiological models with spatial correlations exhibit two 
types of oscillatory behaviour in the endemic phase. In a large parameter range, the oscillations 
are due to resonant amplification of stochastic fiuctuations, a general mechanism first reported for 
predator-prey dynamics. In a narrow range of parameters that includes many infectious diseases 
which confer long lasting immunity the oscillations persist for infinite populations. This efi'ect is 
apparent in simulations of the stochastic process in systems of variable size, and can be understood 
from the phase diagram of the deterministic pair approximation equations. The two mechanisms 
combined play a central role in explaining the ubiquity of oscillatory behaviour in real data and in 
simulation results of epidemic and other related models. 
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Cycles are a very striking behaviour of prey-predator 
systems also seen in a variety of other host-enemy sys- 
tems — a case in point is the pattern of recurrent epi- 
demics of many endemic infectious diseases [ij . The con- 
troversy in the Uterature over the driving mechanisms 
of the pervasive noisy oscillations observed in these sys- 
tems has been going on for long because the sim- 
plest deterministic models predict damped, instead of 
sustained, oscillations. One of the aspects of this contro- 
versy is whether these mechanisms are mainly external or 
intrinsic, and the effects of seasonal forcing termsj^, Q 
and of higher order non-linear interaction terms [5| have 
been explored in the framework of a purely determinis- 
tic description of well-mixed, infinite populations. These 
more elaborate models exhibit oscillatory steady states 
in certain parameter ranges, and have led to successful 
modelling when external periodic forcing is of paramount 
importance @, but they fail to explain the widespread 
non-seasonal recurrent outbreaks found, for instance, in 
childhood infectious diseases |7,]. 

During the last decade, important contributions have 
come from studies that highlight the inherently stochastic 
nature of population dynamics and the interaction pat- 
terns of the population as irnportant endogenous factors 
of recurrence or periodicity [8| . A general mechanism of 
resonant amplification of demographic stochasticity has 
been proposed to describe the cycling behaviour of prey- 
predator systems [9] and applied recently to recurrent 
epidemics of childhood infectious diseases The role 
of demographic stochasticity modelled as additive Gaus- 
sian white noise of arbitrary amplitude in sustaining in- 
cidence oscillations had long been acknowledged in the 
literature [llj . The novelty in [91] and [10|] was that of pro- 
viding an analytical description of demographic stochas- 
ticity as an internal noise term whose amplitude is de- 
termined by the parameters and the size of the system 
using a method originally proposed by van Kampen [l^] ■ 

Our goal is to extend this approach by relaxing the 
homogeneous mixing assumption to include an implicit 



representation of spatial dependence. We show that the 
inclusion of correlations at the level of pairs leads to dif- 
ferent quantitative and qualitative behaviours in a re- 
gion of parameters that corresponds to infectious diseases 
which confer long lasting immunity. Our motivation was 
twofold. On one hand, the homogeneous mixing assump- 
tion is known to give poor results for lattice or network 
structured population [3], 14 1. On the other hand, sys- 



tematic simulations of infection on small- world networks 
have shown that the resonant amplification of stochastic 
fluctuations is significantly enhanced in the presence of 
spatial correlations [15|. Therefore, apart from stochas- 
ticity, the correlations due to the contact structure are 
another key ingredient to understand the patterns of re- 
current epidemics. One of the main difficulties in in- 
cluding this ingredient is that the relevant network of 
contacts for disease propagation is not well known 14|. 
In this paper we shall consider a stochastic Susceptible- 
Infective-Recovered-Susceptible (SIRS) epidemic model 
that leads to the ordinary pair approximation (PA) equa- 
tions of [isj in the thermodynamic limit as the simplest 
representation of the spatial correlations on an arbitrary 
network of fixed coordination number k. The power spec- 
trum of the fiuctuations around the steady state can be 
computed following the approach of [9] and [12]. The 
combined effect of stochasticity and spatial correlations 
has been much studied through simulations, but this is 
an analytical treatment of a model that includes both 
these ingredients. 

Consider then a closed population of size at a given 
time i, consisting of ni individuals of type S, n2 indi- 
viduals of type /, and (A^ — ni — 712) individuals of type 
R, modelled as network of fixed coordination number k. 
Recovered individuals lose immunity at rate 7, infected 
individuals recover at rate 6, and infection of the suscep- 
tible node in a susceptible-infected link occurs at rate A. 
Let ria (respectively, and 715) denote the number of 
links between nodes of type S and / (respectively, S and 
R and R and /). In the infinite population limit, with 
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the assumptions of spatial homogeneity and uncorrelated 
pairs, the system is described by the deterministic equa- 
tions of the standard or uncorrelated PA as follows |13l |: 
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In the above equations the variables stand for the limit 
values of the node and pair densities pi = ni/N, p2 = 
n2/N, p3 = n3/{kN), pi = n^j (kN), and p5 = nb/{kN) 
as N 00. As expected, neglecting the pair correlations 
and setting p3 — piP2 in the first two equations leads to 
the classic equations of the randomly mixed or mean-field 
approximation (MFA) SIRS model. 
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The phase diagrams of the two models are plotted in Fig. 
[T] [solid lines for Eq. ([1]) and dashed line for Eq. ^ , both 
with fc = 4]. We have set the time scale so that ^ = 1. 
The critical line separating a susceptible-absorbing phase 
from an active phase where a stable steady state exists 
with nonzero infective density is given by A^^^ = j 
(dashed black hne) for the MFA, and by AP^(7) = 3^ 
[solid orange (gray) line] for the PA. In addition, in the 
active phase of the PA we find for small values of 7 a new 
phase boundary [solid blue (black) line] that corresponds 
to a Hopf bifurcation and seems to have been missed in 
previous studies of this model [13] ■ This boundary sep- 
arates the active phase with constant densities (region 

I) from an active phase with oscillatory behavior (region 

II) . The maximum of the curve is situated at A ~ 2.5, 
7 ~ 0.03, which means that the PA model predicts sus- 
tained oscillations in the thermodynamic limit when loss 
of immunity is much slower than recovery from infection. 
According to published data for childhood infections in 
the pre- vaccination period Q , taking the average immu- 
nity waning time to be of the order of the length of the 
elementary education cycles at that time (10 years for 
the data points in Fig. [T]) many of the estimated param- 
eter values for these diseases fall into oscillatory region 
II, and the others are in region I close to the boundary. 
Different data points for the same disease correspond to 
estimates for A based on different data records. Although 
small enough to be missed in a coarse grained numeri- 
cal study, the oscillatory phase is large in the admissible 
parameter region of an important class of diseases. A 
systematic study of the dependence of this oscillatory 
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FIG. 1: (Color online) Phase diagram in the (A, 7) plane for 
the MFA and the PA deterministic models and parameter val- 
ues for measles (A), chicken pox (o), rubella (□) and pertussis 
(o) from data sources for the pre-vaccination period. The blue 
stars are the parameter values used in Fig. O 



phase on the parameter k and of its relevance to under- 
stand the behaviour of simulations on networks will be 
reported elsewhere [I6)] . Preliminary results indicate that 
the oscillatory phase persists in the range 2 < fc < 6, and 
that it gets thinner as k increases. There are indications 
that this oscillatory phase is robust also with respect to 
variations of the underlying model 17 1. 

Let us now study the combined effect of correlations 
and demographic stochasticity in region I by taking N 
large but finite. In the stochastic version of the MFA 
SIRS model the state of the system is defined by ni and 
n2 which change according to the transition rates as 
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associated to the processes of immunity waning, recov- 
ery and infection. Here 7^"i+fci."2-i-fc2 denotes the transi- 
tion rate from state (^1,712) to state (ni -I- ki,n2 + k2), 
ki e {-1,0,1}, where i = 1,2. As in the power 
spectrum of the normalized fluctuations (PSNF) around 
the active steady state of system ^ can be computed 
approximately from the next-to-leading-order terms of 
van Kampen's system size expansion of the correspond- 
ing master equation. Setting ni{t) = Npi{t) + \/Nxi{t) 
and 712 (i) = Np2{t) + ^/Nx2{t), the equations of motion 
for the average densities ([2]) are given by the leading- 
order terms of the expansion. The next-to-leading-order 
terms yield a linear Fokker-Planck equation for the prob- 
ability distribution function n(a;i (t) , X2 (t) , t) . The equiv- 
alent Langevin equation for the normalized fluctuations 
is Xi{t) — J2j=i Jij^ji^) +Li{t), where J is the Jacobian 
of Eq. ^ at the endemic equilibrium and L, (t) are Gaus- 
sian white noise terms whose amplitudes are given by the 
expansion. Taking the Fourier transform we obtain for 
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the PSNF, 

P,M ^ (|i,MP) =J2m-\u;)B,,MI^\~u;) , (4) 

where Mik{Lo) = icuSik — Jik and {Li(uj)Lj{uj')) — 
BijSiuj + uj'). For fc = 4 and (5=1 this expression be- 
comes 
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where D and T are the determinant and the trace of J 



and Bii = B 
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susceptible and the infected PSNFs, respectively. 

In a stochastic version of the PA SIRS model the 
state of the system is defined by the integers n^, where 
i = 1,...,5, and recovery, loss of immunity, and infec- 
tion induce different transitions according to the pairs 
or triplets involved in the process. The simplest set of 
transitions and transition rates compatible with Eq. ([l]) 
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This is a coarse grained description where the effect of 
the change in state of a given node on the k pairs that it 
forms is averaged over each pair type. For instance, the 
event of loss of immunity occurs at a rate 777 a, where 
ur is the number of recovered nodes, and changes the k 
pairs formed by the node that switches from recovered 
to susceptible. On average, each pair type will change 
by k units at a rate proportional to its density, accord- 



ing to the equation 777^? 



+ 



kufi kn 

where nRR is the number of pairs of recovered nearest 
neighbours. Taking this level of description and using 
fc7ij^ = 774 -|- 775 -f 277flfl and 771 + 772 + 77^ = A'^, we obtain 
the first three equations of Eq. ([7]) for the rates of the 



three different pair events associated with loss of immu- 
nity. A full microscopic description would require con- 
sidering separately all possible five-node configurations 
for the central node that switches from R to S and its 
four nearest neighbours. We have checked that the de- 
tailed stochastic model involving 40 different transitions 
for fc = 4 gives essentially the same results [l^ as the 
coarse grained model d?]) that we consider here. 

For the fluctuations of the pair densities we set 773 (t) = 
Nkpsit) + VNkxsit), ni{t) = Nkpi{t) + ^Nkx^it), and 
115(1) = Nkp^lt) + y/Nkx^it). The leading order terms 
of van Kampen's system size expansion of the master 
equation associated to ([7]) yield the deterministic PA Eqs. 
([1]). An approximate analytical expression for the PSNF 
can be obtained as before from the next-to-leading-order 
terms. Formula @ is still valid taking now J as the 
Jacobian of Eq. ([1]) at the endemic equilibrium and the 
noise cross correlation matrix B computed directly from 
the expansion. 




FIG. 2: (Color online) (a) Analytical and numerical PSNFs 
of the infectives for model ^ with 7 = 0.1 and A = 2.5; (b) 
the same for model ((Tjl; (c) a similar plot for model ((T)) with 
7 = 0.034 and A = 2.5, notice the lin-log scale; (d) location of 
the parameter values chosen for (a) and (b) (circle) and for (c) 
(square); and (e) and (f) plots of the peak amplitude of the 
PSNF of the PA model as a function of N for the parameter 
values chosen for (b) and (c). 

In Fig. [2] the approximate PSNFs given by Eq. (|4]) 
are plotted (black lines) and compared with the numer- 
ical power spectra of stochastic simulations for N = 10^ 
[green (gray) lines] for models ((3|) and (O [Figs. [2ta) 
and mb)]. For this system size, there is almost per- 
fect agreement between the analytical approximate ex- 
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pressions and the results of the simulations across the 
whole region I. The plots show that for the same param- 
eter values the fluctuations are larger and more coherent 
for model ([7]), in agreement with the results of simula- 
tions reported in [la | for small-world networks on a lat- 
tice and variable small-world parameter. This effect be- 
comes much more pronounced as the boundary between 
regions I and II, where the analytical PSNF of model 
dH) diverges, is approached. Close to this boundary [see 
Fig- M,c)] , there is a significant discrepancy between the 
analytical (black line) and the numerical [green (gray) 
line] PSNFs associated with the appearance of secondary 
peaks at multiples of the main peak frequency. This is 
a precursor of the oscillatory phase, and the breakdown 
of van Kampen's approximation for this system size may 
be understood as an effect of the loss of stability of the 
endemic equilibrium close to the boundary. Relaxation 
towards equilibrium becomes slow compared with the pe- 
riod of the damped oscillations, and a significant part of 
the power spectrum energy shows up in the secondary 
harmonics. For these parameter values, van Kampen's 
expansion becomes a good approximation only for larger 
system sizes. Also shown in Fig. [2je) [respectively, Fig. 
[2jf)] is the scaling with system size of the peak ampli- 
tude of the infectives PSNF of the PA model [pink (black) 
dots] for the parameter values considered in (b) [respec- 
tively, (c)] and the peak amplitude (dashed blue line) of 
the approximated PSNF given by Eq. (|4]). Away from 
the phase boundary of the oscillatory phase we find that 
the simulations exhibit the amplitude and scaling pre- 
dicted by Eq. ([4]) down to system sizes of ~ 10^. By con- 
trast, close to the phase boundary the match is reached 
only for system sizes larger than 5 x 10®. 
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FIG. 3: (Color online) The steady-state infective density 
given by the PA deterministic model (dashed blue lines) , and 
by simulations of the PA and the MFA stochastic models [solid 
black and green (gray) lines, respectively] in regions II and I 
for N = 10^. Parameters are (a) A = 7.5, 7 = 0.01 and (b) 
A = 9, 7 = 0.01. 

Examples of typical time series predicted by the PA 
model in the parameter region of childhood infectious 
diseases are shown in Fig. [3] [dashed blue lines for the de- 
terministic model ^ and solid black lines for simulations 
of the stochastic model The results of simulations 

of the MFA stochastic model ^ for the same parameter 
values are also shown for comparison [solid green (gray) 



lines]. Fig. illustrates the regular high-amplitude 

oscillations of region II. All over this region, simulations 
of the stochastic model ^ reproduce the behaviour of 
the solutions of Eq. ([T]) with added amplitude fluctua- 
tions. The only limitation to observe these oscillations 
in finite systems is that N has to be taken large enough 
for the deep interepidemic troughs to be spanned with- 
out stochastic extinction of the disease. In region I [Fig. 
^h)] there are no oscillations in the thermodynamic limit 
but, in contrast to the stochastic MFA model, the reso- 
nant fluctuations in the PA model are large and coherent 
enough to provide a distinct cycling pattern, which is 
partially described by van Kampen's expansion ([4|. 

In conclusion, we have considered a stochastic version 
of the basic model of infection dynamics including a rep- 
resentation of the spatial correlations of an interaction 
network through the standard PA. We have shown that 
in general the resonant amplification and the coherence 
of stochastic fluctuations are much enhanced with re- 
spect to the MFA model. This quantitative difference 
becomes qualitative in a region of parameter space that 
corresponds to diseases for which immunity waning is 
much slower than recovery. In this region the nonlineari- 
ties of the model and demographic stochasticity give rise 
either to oscillations that persist in the thermodynamic 
limit or to high amplitude, coherent resonant fluctua- 
tions, providing realistic patterns of recurrent epidemics. 

These results are relevant for other population dynam- 
ics models in the slow driving regime that corresponds to 
small 7 in our model, suggesting that in systems of mod- 
erate size intrinsic stochasticity together with the sim- 
plest representation of spatial correlations may be enough 
to produce distinct oscillatory patterns. This favours the 
view that, for a large class of systems, noisy oscillations in 
population dynamics data may be intrinsic, rather than 
externally driven. 
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